---
title: "Latta et al 2016 Costa Rica: Multilevel Regression to determine importance of species traits"
author: "brouwern@gmail.com"
date: "December 22, 2016"
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE, 
                      cache = TRUE)
```


# Intro

This script runs multilevel/hiearchical models with boths months of data (Aug and Jan.)  It calls a number of functions contained in scripts to 

* run a series of nested models, 
* test for sig. of predictors, set up ANOVA tables, and 
* process the output of the models.


Note: The workflow doesn't work perfectly for models that don't 

* have categorical predictors for 2 levels (ie 3-level guild model), 
* a continous predictor (habitat breadth) or 
* data for both months (neotrop migration)

Workflow produces two major outputs

* _model_list.RData = list of the models
* _output_list.RData = ANOVAs of nested models, model list, and metadata; used for final paper



# Preliminaires

## Load Costa Rica data

Set the working directory and load the data.

```{r, message=FALSE, warning=FALSE}
#setwd
my.dir <-  "C:/Users/lisanjie2/Dropbox/Latta_et_al_2017_PeerJ_FINAL"

setwd(my.dir)

#load data
CR.long <- read.csv("./DATA/Latta_et_al_mistnet_long.csv")

```


## Load packages

```{r, message=FALSE, warning=FALSE}
library(lme4) #for multilevel modeling
library(arm)  #for display() command
library(ggplot2)
library(multcomp)
```



## Source my functions for processing models

* There are a number of functions to run the models and process the output.  
* For each species trait (habitat type, diet, etc) a series of models are run (null, time, time*trait) and the output stored in a list.  
* The list is then processed to extract p-valeus, regression coefficients etc. 
* See Latta_et_al_2017_PeerJ_README.Rmd for more information

```{r}
#Main function that runs the models
source("./FUNCTIONS/FX_run_multi_level_model.R")


#Wrapper for helper functions
source("./FUNCTIONS/FX_multlevel_model_wrapper_for_helper_functions.R")


#Misc helper functions for processing model list
#Includes a wrapper that runs all of the functions on the list
source("./FUNCTIONS/FX_multlevel_model_helper_functions.R")

#Function to Extract df, ch2, p from anova() object of two nested modeld
source("./FUNCTIONS/FX_anova_extract.R")

#Round off numeric columns in a dataframe
source("./FUNCTIONS/FX_round_df_columns.R")


```


## Subset of Costa Rica data

* Both months (AUg and Jan)
* Observations seen at least three times (ie "yrs.not.obs" < 7), where yrs.not.obs is the total number of years a species was not observed.

```{r}
i.BOTHMONTHS.obs3.or.more <- with(CR.long, 
                             which(yrs.not.obs < 7))
length(i.BOTHMONTHS.obs3.or.more)#1240
```


Working dataframe to use

```{r}

df <- CR.long[i.BOTHMONTHS.obs3.or.more,]

```



## Run Costa Rica Trait Models


### Run primary habitat model
Model primary habitat: 2ndary forest or primary forest

```{r hab.prime.mod.list, cache=T}

#run primary habitat
## run.multlev.model() calls glmer()
## multiple times
## throws warnings but they
## are not consequential
hab.prime.mod.list <- run.multlev.model(df = df,trait.name = "hab.prime")
  
#save model list
save(hab.prime.mod.list, 
     file = "./MODEL_OUTPUT/hab_prime_model_list.RData")

#process model list
hab.prim.output.list <- process.ML.model.wrapper(model.list = hab.prime.mod.list,
                                                trait.column.name = "hab.prime",
                                                trait.label = "Primary Habitat")  

#save processed list
save(hab.prim.output.list, 
     file = "./MODEL_OUTPUT/hab_prime_output_list.RData")

#look at the ANOVA table
hab.prim.output.list$anova.table
```





### Run disturbance sensitivity model

```{r dist.sense.mod.list, cache = T}
summary(factor(df$sense2))

#disturbance
dist.sense.mod.list <- run.multlev.model(df = df,
                                         trait.name = "sense2")
  
save(dist.sense.mod.list, file = "./MODEL_OUTPUT/dist_sense_model_list.RData")

dist.sense.output.list <- process.ML.model.wrapper(model.list = dist.sense.mod.list,
                                                trait.column.name = "sense2",
                                                trait.label = "Dist. sensivity")  

save(dist.sense.output.list, file = "./MODEL_OUTPUT/dist_sense_output_list.RData")

dist.sense.output.list$anova.table
```





### Run Conservation priority (2 levels) model


```{r cons.priority.mod.list, cache = T}
summary(df$cons.priority)
df$cons.priority <- factor(df$cons.priority)
#disturbance
cons.priority.mod.list <- run.multlev.model(df = df,
                                         trait.name = "cons.priority")
  
save(cons.priority.mod.list, file = "./MODEL_OUTPUT/cons_prior_model_list.RData")

cons.priority.output.list <- process.ML.model.wrapper(model.list = cons.priority.mod.list,
                                                trait.column.name = "cons.priority",
                                                trait.label = "Cons. priority")  

save(cons.priority.output.list, file = "./MODEL_OUTPUT/cons_priority_output_list.RData")

cons.priority.output.list$anova.table
```






### Run Elevational migration  model


```{r elev.mig.mod.list, cache = T}
summary(factor(df$mig.elevYN))
df$mig.elevYN <- factor(df$mig.elevYN)

#disturbance
elev.mig.mod.list <- run.multlev.model(df = df,
                                         trait.name = "mig.elevYN")
  
save(elev.mig.mod.list, file = "./MODEL_OUTPUT/elev_mig_model_list.RData")

elev.mig.output.list <- process.ML.model.wrapper(model.list = elev.mig.mod.list,
                                                trait.column.name = "mig.elevYN",
                                                trait.label = "Elev. mig")  

save(elev.mig.output.list, file = "./MODEL_OUTPUT/elev_mig_output_list.RData")

elev.mig.output.list$anova.table
```






### Run Canopy obligates  model


```{r canopy.ob.YN.mod.list, cache = T}
summary(factor(df$canopy.ob.YN))
df$canopy.ob.YN <- factor(df$canopy.ob.YN)

#canopy.ob.YN
canopy.ob.YN.mod.list <- run.multlev.model(df = df,
                                         trait.name = "canopy.ob.YN")
  
save(canopy.ob.YN.mod.list, file = "./MODEL_OUTPUT/canopy_obYN_list.RData")

canopy.ob.YN.output.list <- process.ML.model.wrapper(model.list = canopy.ob.YN.mod.list,
                                                trait.column.name = "canopy.ob.YN",
                                                trait.label = "Canopy obligate")  

save(canopy.ob.YN.output.list, file = "./MODEL_OUTPUT/canopy_ob_output_list.RData")

canopy.ob.YN.output.list$anova.table
```






### Run canopy use  model


```{r elev.mig.mod.list, cache = T}
summary(factor(df$canopy.use.YN))
df$canopy.use.YN <- factor(df$canopy.use.YN)

#canopy.ob.YN
canopy.use.YN.mod.list <- run.multlev.model(df = df,
                                         trait.name = "canopy.use.YN")
  
save(canopy.use.YN.mod.list, file = "./MODEL_OUTPUT/canopy_useYN_list.RData")

canopy.use.YN.output.list <- process.ML.model.wrapper(model.list = canopy.use.YN.mod.list,
                                                trait.column.name = "canopy.use.YN",
                                                trait.label = "Canopy use")  

save(canopy.use.YN.output.list, 
     file = "./MODEL_OUTPUT/canopy_use_output_list.RData")

canopy.use.YN.output.list$anova.table
```



### Run foraging guild model: 2 levels (specialist vs. omnivore)


```{r elev.mig.mod.list, cache = T}
df$forage.guilde.2lev <- ifelse(df$forage.guilde2 == "F/N" | 
                              df$forage.guilde2 == "I", "specialist","omnivore") 
summary(factor(df$forage.guilde.2lev))

#canopy.ob.YN
forage.guilde.2lev.mod.list <- run.multlev.model(df = df,
                                         trait.name = "forage.guilde.2lev")
  
save(forage.guilde.2lev.mod.list, file = "./MODEL_OUTPUT/forage_guilde2lev_list.RData")

forage.guilde.2lev.output.list <- process.ML.model.wrapper(model.list = forage.guilde.2lev.mod.list,
                                                trait.column.name = "forage.guilde.2lev",
                                                trait.label = "Forage guild")  

save(forage.guilde.2lev.output.list, 
     file = "./MODEL_OUTPUT/forage_guilde2lev_output_list.RData")

forage.guilde.2lev.output.list$anova.table
```













### Run foraging guild model: 2 levels (specialist vs. omnivore)
foraging guilde: 3 levels (forage.guilde2; 2 = 2nd version of factor, NOT 2 levels)

```{r elev.mig.mod.list, cache = T}
summary(factor(df$forage.guilde2))

#canopy.ob.YN
forage.guilde.3lev.mod.list <- run.multlev.model(df = df,
                                         trait.name = "forage.guilde2")
  
save(forage.guilde.3lev.mod.list, file = "./MODEL_OUTPUT/forage_guilde3lev_list.RData")

forage.guilde.3lev.output.list <- process.ML.model.wrapper(model.list = forage.guilde.3lev.mod.list,
                                                trait.column.name = "forage.guilde2",
                                                trait.label = "Forage guild (3 groups)")  

save(forage.guilde.3lev.output.list, 
     file = "./MODEL_OUTPUT/forage_guilde3lev_output_list.RData")

forage.guilde.3lev.output.list$anova.table
```









### Run habitat breadth model: continous numeric variable
foraging guilde: 3 levels (forage.guilde2; 2 = 2nd version of factor, NOT 2 levels)

```{r hab.breadth, cache = T}
#Can take on 7 different values
summary(factor(df$hab.breadth))


hab.breadth.mod.list <- run.multlev.model(df = df,
                                         trait.name = "hab.breadth")
  
save(hab.breadth.mod.list, file = "./MODEL_OUTPUT/hab.breadth_list.RData")

hab.breadth.output.list <- process.ML.model.wrapper(model.list = hab.breadth.mod.list,
                                                trait.column.name = "hab.breadth",
                                                trait.label = "Habitat breadth")  

save(hab.breadth.output.list, 
     file = "./MODEL_OUTPUT/habitat_breadt_output_list.RData")

hab.breadth.output.list$anova.table
```








### Migration status modeling

Lat. / neotrop migrants only occur in January and don't fit into the normal workflow.

### Subset January only

Lat. / neotrop migrants only occur in January
```{r}
i.obs3.or.more.JANUARY.ONLY <- with(CR.long, 
                           which(yrs.not.obs < 7 & month == "Jan"))
length(i.obs3.or.more.JANUARY.ONLY)#560
```

 
Rename dataframe 
```{r}
df.Jan.only <- CR.long[i.obs3.or.more.JANUARY.ONLY,]
```


Set residents as baseline
```{r}
df.Jan.only$status <- factor(df.Jan.only$status, 
                             levels = c("PR","M"))
```


Load function to run models 

This is a modified function that doesn't run mdoels w/ month effets
```{r, cache = TRUE}
source('./FUNCTIONS/FX_run_multi_level_model_NO_MONTH_EFFECT.R')

lat.mig.mod.list <- run.multlev.model.NO.MONTH(df = df.Jan.only, 
                           trait.name = "status")

save(lat.mig.mod.list, file = "./MODEL_OUTPUT/latmig_mod_list.RData")
```



Get output from model
```{r}
str(lat.mig.mod.list,1)

#Summary of best model
summary(lat.mig.mod.list$m7a.best )


#ANOVAs
anova(lat.mig.mod.list$m6.drp.yrXtrt,
      lat.mig.mod.list$m5a.drop.trt)


anova(lat.mig.mod.list$m6.drp.yrXtrt,
      lat.mig.mod.list$m5c.drop.yr)

#ch2 = 0.17
#p =  = .7272


#best (full model)
anova(lat.mig.mod.list$m7a.best,
      lat.mig.mod.list$m6.drp.yrXtrt)

#df 8,9 chi2 2.57, p 0.11

```





## Use multcomp to combine parameters of model

A function to combine individual parameter estimates to get overal tends (ie, combine year + year*slope terms).  Needs to be modified to from previous functions to work with data for lat. migrants b/c they only occur in January.
```{r}
calc.realized.slopes2 <- function(model.list,trait){
library(multcomp)
working.model <- model.list$m7a.best
contrast.mat <-  t(as.matrix(c(0, 0, 0, 1, 1)))

contrast.mat <-  rbind(c(0, 0, 1, 0),
                       c(0, 0, 1, 1))


levs <- levels(working.model@frame$trait)

rownames(contrast.mat) <- levs

summary(working.model)$coefficients

glht.out <- glht(working.model,
     linfct = contrast.mat,
     alternative = c("two.sided"))

glht.out2 <- summary(glht.out,test = adjusted(type = "none"))



#trait.names <- paste(working.list$trait,levs,sep = ".")
trait.names <- trait
df.out <- data.frame(Trait = trait.names,
           Effect.size = glht.out2$test$coefficients,
           SE = glht.out2$test$sigma,
           p = glht.out2$test$pvalues)

df.out$sig <- ifelse(df.out$p > 0.2, "","-")
df.out$sig <- ifelse(df.out$p < 0.10, ".",df.out$sig)
df.out$sig <- ifelse(df.out$p < 0.05, "*",df.out$sig)
return(df.out)
 
}
  

```


Run function to get realized slopes
```{r}

lat.mig.mod.list$m7a.best

lmi <- calc.realized.slopes2(model.list=lat.mig.mod.list, 
                             trait = "Lat. mig")

lmi

#    Trait Effect.size         SE         p sig
# PR  test  0.01977736 0.02003482 0.3235696    
# M   test -0.04833039 0.03831233 0.2071346 
```


Save lat mig object for later use

```{r}
save(lmi,file = "./MODEL_OUTPUT/model_output/lat_mig_multcomp_out.RData")
```

